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Crop models are widely used for the modeling and prediction of crop yields, as decision support tools, 
and to develop research questions. Though typically constructed as a set of dynamical equations, crop 
models are not often analyzed from a specifically dynamical systems point of view, despite its potential to 
elucidate the roles of feedbacks and internal and external forcings on system stability and the optimization 
of control protocols (e.g., irrigation and fertilization). Here we develop a minimal dynamical system, based 
in part on the widely known AquaCrop model, consisting of a set of ordinary differential equations (ODE’s) 
describing the evolution of canopy cover, soil moisture, and soil nitrogen. These state variables are coupled 
through canopy growth and senescence, the evapotranspiration and percolation of soil moisture, and the 
uptake and leaching of soil nitrogen. The system is driven by random hydroclimatic forcing. Important 
crop model responses, such as biomass and yield, are calculated, and optimal yield and profitability 
under differing climate scenarios, irrigation strategies, and fertilization strategies are examined within 
the developed framework. The results highlight the need to maintain the system at or above resource 
limitation thresholds to achieve optimality and the role of system variability in determining management 
strategies. 

© 2017 Elsevier B.V. All rights reserved. 


1. Introduction 

As tools to forecast or backcast crop yields, improve man¬ 
agement strategies, and better understand the physical processes 
underlying crop production, crop models are important tools from 
both a research and an engineering viewpoint (Wallach et al., 2006; 
Steduto et al., 2009). The model outputs, structure, parameteriza¬ 
tion, and data assimilation are all active areas of crop modeling 
research. Because different users have different goals, several types 
of crop models have been proposed, which can be categorized in 
a number of ways. One of the most basic distinctions is between 
dynamic crop models, which are comprised of a set of differen¬ 
tial equations, which are then integrated in time to simulate the 
crop responses of interest at each time point (often daily), and crop 
response models, which, though they may be built on dynamic 
models, relate crop responses directly to inputs (Thornley and 
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Johnson, 1990; Wallach et al., 2006). Most crop models have as 
their main state variables above-ground biomass, leaf area index 
(LAI), harvestable yield, and water and nitrogen balances, though 
the choice and precise number of state variables varies (Wallach 
et al., 2006). Virtually all crop models are process-based, but nec¬ 
essarily involve empirical components, and are of varying levels of 
complexity, depending on the particular goals of the model and on 
the availability of input data. Some are specific to certain crops or 
groups of crops, such as CERES (Ritchie et al., 1998) and AZODYN 
(Jeuffroy and Recous, 1999), while others are more generic, such 
as CROPGRO (Boote et al., 1998), CROPSYST (Stockle et al., 2003), 
STICS (Brisson et al., 2003), and some focus on particular regions 
(e.g., INFOCROP (Aggarwal et al., 2006) for tropical regions). Also 
in the category of generic models, but with a more parsimonious 
framework, is AquaCrop (Steduto et al., 2009). Despite the abun¬ 
dance of crop models which have dynamical systems at their core, 
they are not often analyzed as dynamical systems per se - that is, 
using the wide array of tools and methods provided by dynami¬ 
cal systems theory to understand the mathematical behavior and 
properties of the models (Strogatz, 2014). There are a number of 
potential reasons for this, such as the difficulty of applying these 
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methods to complex models and the aims of modelers, which may 
be focused toward other goals. 

Although they tend to be considerably more complex and 
serve different purposes, crop models share many features and 
describe many of the same processes as do minimal ecohy- 
drological models. The use of such models, which are typically 
formulated as dynamical systems, has provided many insights into 
soil moisture dynamics, plant-water interactions, and nutrient 
cycling (Rodriguez-Iturbe et al., 1999; Porporato et al., 2002,2003; 
Rodriguez-Iturbe and Porporato, 2004). Some features of this type 
of ecohydrological model, such as the parsimonious representa¬ 
tion of processes and stochastic and dynamic coupling between 
state variables, are well-suited to study the feedbacks, nonlineari¬ 
ties, and effect of random hydroclimatic forcing on agroecosystems 
(Porporato et al., 2015). Indeed, the underlying assumptions of 
many dynamic ecohydrological models are better met in agroe¬ 
cosystems than in the natural ecosystems where they are normally 
applied. Such assumptions include homogenous soil depth and 
plant spacing, as well as good drainage, which describe well an 
agricultural field with tillage, uniform crop spacing, and tile drains. 

Various studies have used a dynamical systems framework to 
examine grass ecosystems (Thornley and Verberne, 1989; Tilman 
and Wedin, 1991), grass growth modulated by competition with 
legumes (Thornley et al., 1995) and grazing (Johnson and Parsons, 
1985), forest ecosystems (Thornley and Canned, 1992), forest 
ecosystems under harvest (Parolari and Porporato, 2016), soil 
salinity and sodicity (Mau and Porporato, 2015), and the cycles 
themselves, including feedbacks and nonlinearities (Porporato 
et al., 2003; Manzoni et al., 2004; Manzoni and Porporato, 2007). 
Studying crop models with dynamical systems theory allows for the 
more ready exploration of many interesting aspects of crop sys¬ 
tems, including their stability with respect to parameter change, 
the feedbacks between water, carbon, and nutrient cycling, the 
optimal conditions for growth, and the impact of external inputs 
such as changes in climate patterns and management choices (i.e. 
fertilization and irrigation). 

With the goal of taking advantage of the tools of dynamical sys¬ 
tems theory, in this work we develop a dynamic crop model which 
captures the main crop fluxes and responses of interest without 
being overly complex. The model has three main variables which 
interact dynamically; the canopy cover, the relative soil moisture, 
and the soil nitrogen. The differential equations which account 
for these components are coupled via the crop growth, nitrogen 
uptake and leaching, and evapotranspiration terms. Biomass and 
yield, which are not considered to interact dynamically with the 
other state variables but rather are determined by them, are also 
included as derived variables of agroecologic interest. The model is 
used to examine the crop response to water and nutrient availabil¬ 
ity and varying climatic conditions in order to examine questions of 
optimal fertilization and irrigation and reduction of nutrient leach¬ 
ing. 

Several aspects of the model are derived from AquaCrop 
(Steduto et al., 2009; Raes et al., 2009; Hsiao et al., 2009), which 
is the existing generic crop model that, in addition to its parsi¬ 
mony, can perhaps most easily be viewed as a dynamical system. 
It is also physically based, validated for a variety of crops, and 
widely known. AquaCrop itself is largely based on earlier FAO pub¬ 
lications, in particular through its use of crop coefficients (Allen 
et al., 1998) and in the relation between crop water uptake and 
yield (Doorenbos and Kassam, 1998). The most notable similarities 
between the model developed here and AquaCrop are that canopy 
cover is used rather than the more typical LAI, that evapotranspi¬ 
ration is represented by crop coefficients, and in the dependence 
of the partitioning of transpiration and evaporation on the canopy 
cover. Some key differences involve the soil moisture balance (the 
model developed here makes use of a single vertically averaged 


soil moisture value rather than a soil column consisting of multiple 
layers, and it uses the same soil moisture stress thresholds through¬ 
out) and the nitrogen balance (a balance of total mineral nitrogen 
in the soil is used here rather than the empirical fertility coefficient 
employed in AquaCrop). 

Here a different viewpoint and set of tools is emphasized for 
studying dynamic crop models, and we also aim to place crop mod¬ 
els in a dynamical systems context and to discuss the application of 
the associated methods to crop models. We hope that this contri¬ 
bution will be of interest to both the crop modeling community and 
to researchers in the area of theoretical ecohydrology as a means 
to explore the response of agroecosystems to uncertain climatic 
conditions and optimal management strategies. 

2. Model components 

In this section a dynamical system is constructed which 
describes the interaction of three main components: canopy cover 
C(t), relative soil moisture S(t), and total nitrogen content in the 
soil N(t). We also consider two related variables, namely the crop 
biomass B(t ) and the crop yield Y(t) (hereafter we drop the t- 
dependence of the state variables). The model is interpreted at the 
daily timescale (no diurnal dynamics are considered) and applied 
over the course of a single growing season. It can be forced by ran¬ 
dom rainfall inputs (Rodriguez-Iturbe and Porporato, 2004), and is 
assumed to apply to an agricultural field which is homogenous in 
terms of soil composition, climatic forcing, and management. 

2.1. Canopy cover dynamics 

We define the canopy cover to be the fraction of ground covered 
by a crop. The benefit of using this alternative to the LAI, which 
was also employed by AquaCrop (Steduto et al., 2009), is that it 
combines multiple attributes of the crop canopy into a single, easily 
measured or estimated variable. The rate of change in canopy cover 
is modeled as a balance between the increase due to canopy growth 
and the decrease due to metabolic limitations and senescence, so 
that 

^=C(C,S,N,t)-M(C,t), (1) 

where G is the canopy growth rate, and M is a term which combines 
the effects of metabolic limitation and senescence. The growth rate 
is assumed to be proportional to the rate of nitrogen uptake, U 
(discussed further in Section 2.3), giving 

G(C,S,N,t) = r G .H(C,S,N,t), (2) 

where r G is the canopy cover increase per amount of nitrogen taken 
up (the value for this and other crop growth parameters can be 
found in Table 1 ). The combined metabolic limitation and mortal¬ 
ity/senescence term is 

M(C, t) = (r M + y(t - tsen)- ©(t - t sen )) ■ C 2 , (3) 

where the first term, r M , is a constant metabolic limitation term, and 
the next term is a time-dependent mortality and senescence term. 
For the latter, a linear function is used which increases with a slope 
of y after the senescence onset time, t sen , at which point the Heavi¬ 
side step function, 0, causes the senescence term to begin to affect 
the equation. This form recalls somewhat the Gompertz-Makeham 
law (Makeham, 1860), which includes an age-independent mortal¬ 
ity term and an age-dependent mortality term, although here the 
constant term is conceptualized as a metabolic limitation term and 
the time-dependent term as a senescence term. For unstressed con¬ 
ditions (sufficiently high 5 and N) prior to t sen , Eq. (1) is the logistic 
growth equation (Murray, 2002), and it includes the approximately 
exponential growth of C in the initial growth stage, the slowing of 
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Table 1 

Model parameters. 


Parameter 

Value 

Units 

Name/description 

Source 

r G 

560 

m 2 /kg N 

Canopy growth per unit N uptake 

Calculated using data from Hsiao et al. (2009) 

r M 

0.2 

1/d 

Canopy decline due to metabolic limitation 

Calculated using data from Hsiao et al. (2009) 

y 

0.005 

1/d 2 

Slope of increase of senescence after t se n 

Calculated using data from Hsiao et al. (2009) 

I<cb 

1.03 

- 

Max. T/ETq 


Kce 

1.1 

- 

Max. E/ETo 

Hsiao et al. (2009) 

tsen 

110 

d 

Days until onset of senescence 

Mean of values in Table 2 of Hsiao et al. (2009) 

tcs 

140 

d 

Length of growing season 

Mean of values in Table 2 of Hsiao et al. (2009) 

W* 

3.37 xl0“ 2 

kg B/m 2 /d 

Normalized daily water productivity 

Hsiao et al. (2009) 

h 

0.5 

kg Y/kg B 

Maximum harvest index 

Hsiao et al. (2009) 

Vc 

0.054 

kg N/m 3 water 

Maximum N concentration taken up 

Derived from model parameters 

D 

5.5 x 110“ 6 

kg/m 2 /d 

N deposition rate 

National Atmospheric Deposition Program (NRSP-3) (2017) 

Ft 

0.0286 

kg N/m 2 

Maximum N uptake 

Bender et al. (2013) 

Py 

0.12 

S/kg 

Corn price per kg of yield 

Lamm et al. (2007) 

Pf 

0.639 

S/kg 

Fertilizer unit price N 

Lamm et al. (2007) 

Pi 

0.0148 

S/m 3 

Irrigation water unit price 

Vico and Porporato (2011b) 

Pl 

0 

$/kg 

Cost of leached N 

Set to 0 in current simulations 

Pfix 

0.109 

S/m 2 

Fixed costs 



growth as a limit is reached, and the negligible growth rate near 
the carrying capacity. This compares well with the data for canopy 
cover presented by Hsiao et al. (2009) (see Section 3.1 ). 


2.2. Soil moisture balance equation 


the model should be interpreted at the daily timescale. The water 
stress coefficient is given by 


US) = 


o s <s w , 

S — S w C c 

— --— Sw < s < s*, 

J * — Jvv 


( 6 ) 


Soil moisture is modeled as a balance between gains from rain¬ 
fall and irrigation and losses to evapotranspiration and leakage 
(Rodriguez-Iturbe and Porporato, 2004; Vico and Porporato, 2010) 

4>Z^= m + 1(S, t) - T(S, C, t) - E(S, C, t ) - L(S), (4) 

where S is the vertically averaged relative soil moisture, 0 is poros¬ 
ity, and Z is a soil depth with homogenous characteristics (Table 2 
contains values for the soil parameters). 0Z is defined as the active 
soil depth (Laio et al., 2001b), the volume per surface area available 
for water storage. In agricultural soils, tilling tends to rearrange soil 
profiles so that the top layer of soil is relatively uniform in compo¬ 
sition and depth. We assume that the root growth (which we do not 
explicitly model) is constricted to Z, and that hydraulic redistribu¬ 
tion over this depth allows water to easily move to areas of lower 
soil moisture, making the vertically-averaged soil moisture a good 
description of the amount of water available for evapotranspiration 
(Guswa et al., 2002). 

R is the rainfall rate. For the purposes of a probabilistic analysis, 
here it is modeled as a marked Poisson process with mean event 
frequency X and exponentially distributed rainfall events depths a 
(Rodriguez-Iturbe and Porporato, 2004). This stochastic component 
allows for the model to include the effect of unpredictable external 
forcing via rainfall, which is especially important in arid and semi- 
arid ecosystems, and for rain-fed agriculture. 

In the case of irrigated agriculture, a term I gives the irrigation 
rate, which may be a function of S and/or t depending on the irri¬ 
gation strategy employed (e.g., stress avoidance or microirrigation) 
(Vico and Porporato, 2010, 201 la, b). 

The transpiration rate T is assumed to be proportional to C and 
is given by 


T(S, C, t) = K s (S) ■ C K cb ETo(t), (5) 

where K s {S) is a water stress coefficient, U is a basal crop coeffi¬ 
cient (essentially the midseason basal crop coefficient of Allen et al. 
(1998)), and ET 0 is the reference evapotranspiration, which is cal¬ 
culated using the Penman-Monteith equation for a reference crop 
(normally grass, but occasionally alfalfa) (Allen et al., 1998). As no 
diurnal variation is considered, ET 0 is a mean daily rate and thus 


^ 1 > S* < S, 

where 5 W is the wilting point and S* is the point of incipient sto- 
matal closure. I< s {S) therefore captures the plant stomatal response 
to soil moisture conditions. As mentioned previously, the plant is 
assumed to be able to easily compensate for areas of low soil mois¬ 
ture in the soil column by drawing more water from areas of high 
soil moisture, making S a good indicator of the amount of water 
available to the plant. However, this assumption is weakened if 
the plant cannot do so because of high root resistance or spatial 
heterogeneities in the soil properties (Guswa et al., 2002). 

The evaporation rate E is assumed to be proportional to (1 - C) 
and is given by 


E(S, C, t) = US) • (1 - C) • K ec • ET 0 {t), 


(7) 


where US) reduces evaporation according to soil moisture and I< ec 
is a baseline evaporation coefficient. A similar dependence of evap¬ 
oration on 1 - C was used by Steduto et al. (2009). The evaporation 
reduction coefficient is given by 


US) 



( 8 ) 


where S h is the hygroscopic point, below which no soil moisture 
losses occur. A diagram of I< s and I< r as a function of 5 is shown in 
the upper panel of Fig. 1 , and the increase of evapotranspiration as a 
whole with increasing S can be seen from top to bottom in the lower 
panel of Fig. 1. Evaporation draws primarily from a thin top layer of 
soil, drawing from lower soil layers only when potential gradients 
drive water from lower soil depths upward. This is often modeled 
using the two stage method for soil evaporation (Ritchie, 1972; 
Brutsaert and Chen, 1995). The dependence of E on the average soil 
moisture value over a depth Z simplifies the actual relationship, but 
it does capture the high rates of evaporation at saturation (S = 1) and 
the trend toward a rate of zero evaporation as S approaches S^. The 
form that is used for I< s is essentially equivalent to the expression 
for transpiration used in Laio et al. (2001b), while the form for I< r 
is quite different from that used for evaporation in the same paper. 
Laio et al. (2001b) considered evaporation and transpiration sep¬ 
arately, with the former being very small due to the presence of 







Table 2 

Climate and soil parameters. 
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Parameter 

Value 

Units 

Name/description 

Source 

a 

1.5 

cm 

Mean rainfall depth 

Sample values 

X 

0.3 

1/d 

Mean rainfall frequency 

Sample values 

ET 0 

5xl0- 3 

m/d 

Reference evapotranspiration 

Approximated from Hsiao et al. (2009) 

s h 

0.14 

- 

Hygroscopic point 

Rodriguez-Iturbe and Porporato (2004) 

Sw 

0.17 

- 

Wilting point 

Rodriguez-Iturbe and Porporato (2004) 

S* 

0.35 

- 

Point of incipient stomatal closure 

Rodriguez-Iturbe and Porporato (2004) 

Sfc 

0.59 

- 

Field capacity 

Rodriguez-Iturbe and Porporato (2004) 

ksat 

0.33 

m/d 

Saturated hydraulic conductivity 

Rodriguez-Iturbe and Porporato (2004) 

d 

13 

- 

Leakage parameter 

Rodriguez-Iturbe and Porporato (2004) 

a 

~1 

- 

Fraction of N dissolved 

Porporato et al. (2003) 


0.43 

- 

Soil porosity 

Rodriguez-Iturbe and Porporato (2004) 

Z 

1.0 m 

d 

Soil depth 

Irmak and Rudnick (2014) 



SH 

xIO ' 3 



C[-] 


Fig. 1 . Top: the water stress coefficient (dashed line) and the evaporation reduction 
coefficient (solid line) as a function of soil moisture S. Bottom: evapotranspiration 
[m/d] as a function of S and C, with values of ET 0 and the soil moisture thresholds as 
in Table 1. 

the plant canopy. However, as we are interested in the crop canopy 
as it develops throughout the growing season (from left to right in 
the lower panel of Fig. 1), the maximum values for T and F must 
be of somewhat similar magnitude to capture the dominance of E 
shortly after planting and that of T later in the growing season (this 
is reflected in the fact that I< cb and I< ce are indeed nearly the same) 
(Kelliher et al., 1995). 

The combined percolation and runoff rate is denoted as Q, and 
as we are considering well-drained agricultural fields, subsurface 
percolation is assumed to dominate compared to overland runoff 
and to be equal to the hydraulic conductivity, i.e., 

Q(S) = k(S) = ksacS d , (9) 


where k is the hydraulic conductivity, k sa t is the saturated hydraulic 
conductivity, and d is an empirically based parameter (Brooks and 
Corey, 1964; Rodriguez-Iturbe and Porporato, 2004). 

2.2.1. Calculation ofS w and S* 

Using data for silty loam (a common agricultural soil) and meth¬ 
ods from Clapp and Hornberger (1978) and Laio et al. (2001a), 
the wilting point S w of was calculated as the soil moisture level 
corresponding to a matric potential of -1.5 MPa. Corn begins to 
suffer water stress when approximately 50% of the total available 
water (which is the water content at field capacity minus that at 
the wilting point) is depleted (Rhoads and Yonts, 2000). There¬ 
fore, we calculate the point of incipient stomatal closure S* as 
5* = (S w +Sf c )/2. For silty loam, S w = 0.35,5*=0.47, and Sf c = 0.59. 

2.3. Soil nitrogen content 

While AquaCrop (Steduto et al., 2009) makes use of an empiri¬ 
cal measure of soil fertility that allows the model to be used even if 
detailed soil nitrogen data are not available, most crop models con¬ 
sider a nitrogen balance (Ritchie et al., 1998; Jeuffroy and Recous, 
1999; Boote et al., 1998; Stockle et al., 2003; Brisson et al., 2003; 
Aggarwal et al., 2006) due to its key role in the growth and devel¬ 
opment of crops. In order to better examine crop growth and yield 
under optimal fertilization and irrigation strategies, a soil nitro¬ 
gen balance is also included here. The evolution of total mineral 
nitrogen in the soil is given by the balance between deposition 
and fertilization as inputs and leaching and plant uptake as outputs 
(Porporato et al., 2003) 

^ = D(t) + F(N, t) - I(S, N) - U(S, N, C, f), (10) 

where N is nitrogen content per unit area of soil, D is the rate of nat¬ 
ural nitrogen addition to the soil, and Fis the fertilization rate. For all 
figures in this paper, the average annual rate of nitrogen deposition 
for a heavily agricultural region has been used as a constant depo¬ 
sition rate D (National Atmospheric Deposition Program (NRSP-3), 
2017). Unless otherwise noted, the fertilization rate Fis considered 
to be the maximum potential uptake of nitrogen F t divided by the 
length of the growing season, t GS . The total mineral nitrogen con¬ 
tent in the soil, rather than the individual nitrate and ammonium 
components, is used because plants are able to take up both forms, 
making the separation of the two unnecessary in the case of this 
model, which aims for a general picture of nitrogen fluxes. 

The leaching term L is proportional to the percolation from the 
hydrologic balance, Q, and the nitrogen concentration as 

(ID 

where a is the fraction of N which is dissolved in the soil moisture 
(a ^ 1 for nitrate, while a < 1 for ammonium). The nitrogen concen- 
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tration in the soil moisture is given by the quantity which is 
denoted by q. 

Plant uptake of nitrogen, U, is given by 


U(S, N, C, t) = f(q) • T(S, C, t), 


( 12 ) 


in which f[q) is a function which limits the nitrogen uptake above 
a certain critical concentration q c , with the form 


m = 


' aN 

aN 

w 

s$z < f]c ’ 

y 

aN 

Be 

s$z ~ r,c ' 


(13) 


The physical reasoning for this limitation is that beyond a cer¬ 
tain point, taking up more nitrogen is not useful for the plant to 
increase its growth rate, and extremely high nitrogen concentra¬ 
tions in plant tissue are toxic to the plant. The above limitation 
is meant to account in a parsimonious way for the plant’s abil¬ 
ity to exclude nitrogen from transpired water (i.e. active uptake 
(Porporato et al., 2003)). It is worth noting that a reduction in S 
can either increase or decrease the N uptake. As long as 5 >5*, a 
reduction in S increases the concentration q, thereby increasing N 
uptake if initially q < r \ c . However, if S drops below 5*, transpiration 
decreases and therefore so does N uptake. 



Fig. 2. Growth of canopy cover in the C-only model, using the parameterization as 
described in the text (black line). Data are from 6 years of maize experiments in 
Davis, CA (Hsiao et al., 2009) (open circles). 


deterministic dynamical system in order to demonstrate some of 
the insights which can be gained from this approach. 


2.4. Crop biomass and yield 


While the dynamics of the model are contained in the equations 
for C, 5, and N, other variables which depend on one of the three 
main variables are also of interest. Specifically, we consider the crop 
biomass B and crop yield Y. 

The accumulation of plant biomass is modeled using the nor¬ 
malized daily water productivity W* (e.g., Steduto et al. (2009)), 
which is typically multiplied by the ratio of transpiration to ref¬ 
erence evaporation to model biomass accumulation. However, in 
place of transpiration we us the nitrogen uptake divided by the 
nitrogen uptake threshold q c , giving 


dB U{S,N,C,t ) 

dt q c ET 0 (t) 


——K s (S)K c ijf(ri)C. 

Be 


(14) 


The use of ^ rather than T allows us to extend the concept of 
water productivity to also consider the effects of nitrogen limita¬ 
tion. When r] > rjci one recovers the biomass growth equation used 
by Steduto et al. (2009) and others, which considered transpiration 
rather than nitrogen uptake for biomass accumulation. 

The biomass and yield are related through a harvest index, h, 
which is the fraction of the biomass which makes up the yield. 
The harvest index is often modeled as an increasing function in 
time which is modulated by various stresses (Steduto et al., 2009; 
Raes et al., 2009). Here we instead utilize a reference value for h 
and assume that stress limitations are sufficiently accounted for 
elsewhere in the crop growth equations, recognizing that this limits 
the validity of the crop yield calculations to the end of the growing 
season. The yield is then 


Y — h-B. 


(15) 


3. Reduced versions of the model 

The complete model is defined by the balance equations for C, S, 
andN(Eqs. (1), (4), and (10)) and their component fluxes, with Eqs. 
(14) and (15) defining additional variables. Two reduced versions 
of the model will now be examined. The first, in which S and N 
are held constant, is compared to canopy cover data to estimate 
parameters for Eq. (1). The second, in which only S is held constant 
(and there is therefore no stochastic forcing), is analyzed as a typical 


3.1. Canopy growth equation and its parameterization 

We begin by examining the simplest version of the model, in 
which S and N are fixed (at 5 > 5* and with N such that q > q c ) but C 
is allowed to vary. In these conditions (and also with ET 0 constant, 
to maintain analytical tractability), Eq. (1) reduces to 

^ = r G K cb ET 0 r, c ■ C - (r M + y(t - t sen ) ■ ®(t - t se „)) C 2 , (16) 

which is simply the logistic model if t< t se n, in which case this equa¬ 
tion can be solved analytically as 


C(t) = 


rcKchETorjcCoe^^y^ 

r c K cb ET 0 ric + C 0 r M (e^^oict _ ’ 


(17) 


which is the logistic equation (Murray, 2002). 

In order to parameterize Eq. (16), it was necessary to use data 
from a growing season in which the crop did not experience water 
or nitrogen stress. Values for r G , tm, and y have therefore been 
obtained by minimizing the RMSE of the model compared to the 
data from Hsiao et al. (2009) for fully irrigated and fertilized condi¬ 
tions. The data come from 6 seasons of experiments spread over 22 
years in Davis, CA. The first three years used slightly different maize 
cultivars, while the last three used the same cultivars, but in order 
to include more data they have been assumed to be similar enough 
to consider together. An approximate mean reference evapotran- 
spiration rate ET 0 and the value for I< cb were also taken from Hsiao 
et al. (2009). These experiments did not report soil N or uptake 
rates, and so q c has been estimated by averaging the cumulative 
N uptake for maize found in Bender et al. (2013) across the grow¬ 
ing season. While it would be preferable to use a more complete 
single dataset for the model parameterization, the emphasis here 
is not on predicting crop growth but on reproducing the general 
crop behavior. Fig. 2 shows the model vs. the data against which 
it was parameterized, demonstrating a good fit particularly prior 
to the time of senescence, t sen . The value for t se n is cultivar-specific 
and for this figure was taken as the average t se n over the 6 years of 
experiments. The value of this parameter and all others discussed 
in this section can be found in Table 1. 
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Fig. 3. Top: timeseries of canopy cover and soil nitrogen for the C-N model, for two 
differing initial conditions of N (solid lines represent C; dashed lines represent N). 
Bottom: C-N vector plot and phase portrait, for two initial conditions. The black dot 
represents the stable fixed point (N 2 , C 2 ), and the solid gray line is the rj c threshold. 
Note that in this figure, the combined fertilization and deposition rate has been 
reduced slightly in order to better show the impact of varying the initial conditions, 
and the simulation has only been performed until t sen because after this point in 
time the trajectories no longer converge towards the same fixed point. 


3.2. N and C system 

An interesting 2-D dynamical system is obtained when C and 
N are free to vary in time, but S is kept constant, which approxi¬ 
mates the conditions in an agricultural field with a microirrigation 
system and constant fertilization and deposition rate, F+D = Fq. The 
top panel of Fig. 3 shows the evolution of the two state variables 
N and C in time, while the bottom panel is a phase space diagram 
which shows sample trajectories in the C-N phase space. It is easy 
to see the development of the state variables for different initial 
conditions, and the effects of parameter changes on the vector field, 
which determines the direction the system moves for a given condi¬ 
tion, can also be examined using this type of diagram. In the bottom 
of Fig. 3, the ri c threshold of Eq. (13) can be seen as the solid gray 
line which separates the two parts of the solution - on the left side, 
r] < T] c , and on the right, p > p c . Optimization will be further discussed 
in a later section, but for now we point out that in order to max¬ 
imize crop growth, the system should be kept on the right side of 
this threshold, as the trajectories of the vector field here point to 
higher values of C and thereby greater rates of crop growth. 

Different solutions exist above and below p c because when 
p > Pc, sufficient nitrogen is available for crop growth, and Eq. (1) is 
decoupled from N. An analytic expression for C(t) can be obtained 
(when t<t sen ) due to this decoupling, which is shown in Eq. (17). 
An exact expression can also be found for N(t), but as it is rather 
involved it is not included here. When p < p c , the crop experiences 


nitrogen stress and Eq. (1) is again coupled to N. Analytical expres¬ 
sions for C(t) and N(t) are unavailable in this case. 


3.2.1. Fixed points and stability, p>p c 

For the simpler case of p > p c and t< t sen , the first fixed point is 
given by 


C* i = 0, 


N* i 


Fq ■ S(f)Z 
a ■ k sa tS d ’ 


while the expressions for the second are 
C *2 = (PK cb ET 0 r,c . 

r M 

^ = Scpz(F 0 -^Cf) = Scf>Z(F 0 -^ b ETW c ) 
a ■ k sat S d a ■ k sa tS d 


( 18 ) 

(19) 


( 20 ) 

( 21 ) 


Recalling Eq. (1) with ^ = 0, we note that the steady-state uptake 
of nitrogen is given by ^(CJ 2 ) = ^K^ET^p^, a quantity which can 
also been seen inside the parentheses in Eq. (20). The first fixed 
point is an unstable node and the second is a stable node (a third 
exists, but it is always negative and thus not physical). The eigen¬ 
values of the first fixed point are 


_ a ' ksatS d 
la - S$Z~ 


( 22 ) 


b — ETop c K c b r G, 


(23) 


while those of the second fixed point are 

} _ _ d ■ k S atS d 
X2a ~ S0Z 

^2 b = -ETop c K c b r G- 


(24) 

(25) 


The first fixed point is always an unstable node, and in the second is 
always a stable node. This is unsurprising, as the standard logistic 
equation, which is contained within the system dynamics when 
p > Pc, likewise has one stable and one unstable node as its fixed 
points. 


3.2.2. Fixed points and stability, p<p c 

If the system were allowed to develop to steady state (t^ oo), 
the explicitly time-dependent part of the mortality M(C, t) term 
would ultimately drive the canopy cover to a value of zero, and the 
soil nitrogen content would approach a value determined by the 
balance between the fertilization/deposition and leaching terms. 
The fixed points are obtained assuming no senescence term (e.g., 
if a perennial crop rather than an annual one were considered). In 
this condition, there are two fixed points. The first has the same 
expressions as Eqs. (18) and (19), while the second is 


C*9 = 




2 ET 0 K t 


cb 


(26) 


- r M Z(f>k sat S d + 1 + S(f)Z\/(r M k sat S d ) 2 + 4F 0 ET 2 0 I^ b r G 

N * 2 = --• ( 27 ) 

2aFTo I< 2 b r G 


The stability can more easily be seen when it is put in terms of 
the fertilization and deposition term, the critical value of which is 
derived from the expression for the eigenvalues and is given by 

F c = 2aS 2d k] at r M {aETK cb - SZ<pr M ) 

• (a 2 ET 2 K 2 b + 4aETSZ<t>K cb r M -S^tfr 2 ,) (28) 

+ETK cb r c (a 2 ET 2 K? b - 6aETSZct>K cb r M + S 2 Z 2 <t> 2 r 2 M ) 2 , 






































86 


N. Pelak et al. / Ecological Modelling 365 (2017) 80-92 



0 20 40 60 80 100 120 140 


t[d] 

(a) 



(b) 



t [d] 



0 20 40 60 80 100 120 140 

t[d] 


(c) 


(d) 


Fig. 4. Time series of (a) canopy cover C, (b) soil moisture S and rainfall R, (c) soil nitrogen N, (d) crop biomass and yield over a growing season. 


so that the fixed point is a spiral when F 0 < F c and a node when 
F 0 >F C , pointing to the possibility of damped oscillations. Oscilla¬ 
tions related to nitrogen cycling were also observed by Thornley 
et al. (1995) (oscillations of LAI and soil nitrogen in a model of 
grass-legume dynamics), Tilman and Wedin (1991) (oscillations 
and possible chaos in interannual dynamics of a perennial grass), 
Manzoni and Porporato (2007) (shifts in stability in a model of sub¬ 
strate carbon and nitrogen dynamics), and Parolari and Porporato 
(2016) (stability shifts in a model of forest carbon and nitrogen 
cycles under harvesting). The interpretation of such oscillations is 
not entirely clear, as it is possible that they are merely artifacts 
of simplified models or the results of overfitting the available data, 
but their presence in such models is intriguing and deserves further 
attention. 

4. Soil moisture dynamics and hydrologic forcing 

The addition of the soil moisture dynamics greatly increases the 
model complexity, especially when the rainfall stochastic forcing is 
considered. This forcing adds considerable interest to the dynamics 
of the model, as it allows us to consider the effect of varying rain¬ 
fall parameters as well as to examine the model in a probabilistic 
sense. While it is possible to obtain some analytic results regarding 
soil moisture probability distributions for statistically steady states 
under stochastic rainfall forcing (see for example Rodriguez-Iturbe 
and Porporato (2004)), the complexity of the crop growth function 
and nitrogen balance employed here make it necessary to proceed 


numerically (though see Schaffer et al. (2015) for a special case in 
which analytical results were obtained for stochastically driven soil 
moisture and plant biomass). 

4. 1. Soil moisture dry-down 

Many important agroecosytems have some form of the Mediter¬ 
ranean climate, in which the rainfall occurs out of phase with the 
growing season. In this case, the soil moisture dynamics occur as 
a deterministic dry-down, with the exception of whatever small 
amounts of precipitation may occur during the growing season. 
Therefore, all other factors being equal, the crop yield of rain- 
fed (i.e., non-irrigated) agriculture in this type of climate depends 
greatly on the initial condition of soil moisture that is available at 
the beginning of the growing season. Of course, this initial supply 
may also be supplemented by irrigation, which is similar to the case 
considered in Section 3.2. 

4.2. Stochastic forcing 

Figs. 4 and 5 show the development of the three main state 
variables and their associated fluxes over the course of a grow¬ 
ing season, with t = 0 taken as the start of the growing season. Note 
that in these and the proceeding figures, a constant rate of nitrogen 
fertilization/deposition was imposed. As compared to the deter¬ 
ministic scenarios discussed in Section 3, the variables in the full 
model show much greater variability, due to the direct dependence 
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Fig. 6. A sample trajectory shown in the 3-dimensional phase space of C, S, and N 
(black line), and projections onto the three planes (gray lines). Dashed lines denote 
the dynamics which occur after t sen , and the dotted lines the S* threshold (S-C and 
S-N planes) and the r} c threshold (S-N plane). 

of the fluxes on the stochastically driven soil moisture balance. Fig. 6 
shows the main dynamic variables in the three-dimensional phase 
space. Observing this sample time series, excursions below the soil 
moisture threshold S* (the dotted lines that are perpendicular to 
the S axis) and below the p c threshold (the diagonal line on the 


S-N plane) can be seen to coincide with reductions in C. The case 
of water stress depends on S only, while the latter case of nitrogen 
stress involves the interaction of S and N because of their joint effect 
on th ef[r]) limitation function. 

4.3. Impact of rainfall regimes on rain-fed agriculture 

The timing and amount of rainfall exerts a strong control on 
crop growth in rain-fed agriculture. We first examine the effect 
of different rainfall regimes when associated parameters (mean 
event depth a and mean frequency A.) are constant throughout a 
growing season, which is a reasonable approximation for growing 
season conditions in many regions of the world. Despite the fact 
that the parameters are constant in time, there remains a strong 
intra-seasonal time dependence, primarily due to the growth of 
canopy cover, C. This is due to both its growth in time and more 
explicitly through the time-dependence of the M(C, t) term. 

This pattern in time can be seen not only in Figs. 4 and 5 but also 
in Fig. 7, which shows the ensemble average over many simulations 
of canopy cover (top left), soil nitrogen (top right), soil moisture 
(bottom left), and nitrogen leaching (bottom right). In each simu¬ 
lation, the mean rainfall rate was kept constant but a and X were 
changed, to demonstrate the interaction of rainfall event frequency 
and event size. Simulations with larger, less frequent events are 
characterized by reduced canopy cover, and higher rates of nitro¬ 
gen leaching. Flowever, there are also slightly higher levels of soil 
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(a) (b) 




(c) (d) 

Fig. 7. A comparison of the mean (a) canopy cover, (b) soil nitrogen, (c), soil moisture, and (d) leaching across different precipitation regimes for A = 0.1 d -1 (light gray), 
A = 0.3 d -1 (gray), A = 0.5 d -1 (black), with a altered to keep a constant mean rainfall rate of 4.5 mm/d for all figures. A typical fertilization treatment for corn has been applied, 
resulting in the observed jump in N. 


nitrogen which remain, as the reduced soil moisture and canopy 
cover led to a low nitrogen uptake and thus higher soil moisture 
levels. The system shown in Fig. 7 undergoes a typical fertilization 
schedule for corn, in which some fraction £ of the total fertiliza¬ 
tion F t is applied at the beginning of the growing season, and the 
remainder is applied after a period r, resulting in the jump which 
can be observed in N in Fig. 7b. Flere, £ = 0.3 and r = 40 d (see Section 
5.1 and Fig. 10 for a discussion of the optimization of £ and r). 

5. Optimal strategies 

Crop models represent an important tool to study the impact of 
different management strategies aimed at maximizing yield, mini¬ 
mizing water and fertilizer use, reducing the leaching of fertilizers, 
and optimizing the timing of irrigation and fertilization treatments 
under hydroclimatic variability (Wallach et al., 2006). Toward this 
goal we develop a first-order objective function, which considers 
the profit from the sale of produce, costs of fertilizer and irrigation, 
‘environmental cost’ of nitrogen leaching, and fixed costs, 

Pnet = Py ■ y{t GS ) - Pf ■ Ftot - Pi ■ hot - Pi ■ Ltot - Pfix, (29) 

where p Y [S/kg Y] is the unit sale price of the crop yield at 
the end of the growing season, Y(tcs). Pf [S/kg N] and p/ [$/m 3 ] 
are the unit prices of fertilizer and irrigation water, respectively, 


while the cumulative fertilization and irrigation are given by F to t = 
f* GS F(N, t)dt and hot = Jq GS KS, t)dt. The unit ‘environmental cost’ 
of leached nitrogen is given by p L [$/kg N], here conceptualized 
as the cost necessary to mitigate these losses or to pay associ¬ 
ated fines, while the cumulative nitrogen leaching is given by 
L tot = Jq GS L(S,N)dt. Finally, p^ x [$/m 2 ] is a fixed cost represent¬ 
ing distribution and energy costs, here estimated following Vico 
and Porporato (2011b). Estimated values for these parameters can 
be found in Table 1. This objective function should be thought of 
as a means to quantify the relative financial impact and impor¬ 
tance of various components of the crop system, rather than as a 
way to obtain firm predictions about the profitability of various 
management strategies. 

Fig. 8 shows several key crop responses under idealized, non¬ 
stochastic conditions. The responses of crop yield Y and the 
objective function P net to different mean soil moisture and soil 
nitrogen conditions are illustrated in Fig. 8a and b. Fig. 8c and 
d shows the cumulative amounts of irrigation and fertilization, 
respectively, that would be needed to keep the S and N at the des¬ 
ignated mean values. The horizontal lines mark the 5* threshold, 
below which point the transpiration begins to decrease, and the 
diagonal line marks the p c concentration threshold. 

Nitrogen use efficiency (NUE) is the ratio of the amount of nitro¬ 
gen which is taken up by the crop to the amount which is applied, 
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Fig. 8. Response to the necessary rates of irrigation and fertilization to keep S and N at the designated constant values of (a) crop yield, (b) profit, (c) cumulative fertilization 
[kg/m 2 ], (d) cumulative irrigation [m], (e) nitrogen use efficiency (NUE), the cumulative nitrogen uptake as a fraction of total fertilization, and (f) irrigation efficiency (IE), the 
cumulative transpiration as a fraction of the total irrigation. The horizontal dashed line represents the S* threshold, while the diagonal dashed line is the rj c limit. All figures 
are for deterministic conditions (i.e., no stochastic forcing in the rainfall). 


and is an important metric by which to judge fertilization strategies. 
Similarly, we can define the irrigation efficiency (IE) as the ratio of 
the irrigation water applied to the amount which is transpired by 
the crop. These two values are shown in Fig. 8e and f, respectively, 
as a function of the mean soil moisture and mean soil nitrogen. 
Each of the panels in Fig. 8 sheds light on a different consideration 


for the optimal use of water and nitrogen resources. However, the 
common thread between them is that in each case, the ‘best’ sce¬ 
nario from the perspective of using water and nitrogen resources in 
the most efficient manner (i.e. maximizing IE and NUE to produce 
the highest possible yield and profit), occurs at the intersection of 
the 5* and r] c lines. At this point neither water nor nitrogen is lim- 
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Fig. 9. Yield (top) and profit (bottom) as a function of the constant fertilization rate. 
The dashed line represents a theoretical maximum while the solid line is the mean of 
many simulations. The inset plots show the numerical probability density functions 
at each point. 


iting, and no extra irrigation or fertilization beyond what is needed 
to keep the system at these 5 and N values is used. However, with 
the addition of the random rainfall in the next section, additional 
concerns such as the robustness of the optimal strategies under 
stochastic forcing must also be considered. 

5.1. Optimization under stochastic rainfall conditions 

Random hydroclimatic forcing adds uncertainty to the expected 
value of the objective function. This is illustrated in Fig. 9, which 
shows the numerical probability distribution functions of yield 
and profit for varying fertilization rates. Note that as the fertil¬ 
ization rate increases, both the mean of the yield and its variance 
increase. The mean increases because the higher fertilization rates 
lead to less likelihood that the crop will experience a shortage of 
nitrogen, while the variance increases because the field of possible 
yields expands - the extra nitrogen raises the maximum possible 
yield, while a low-rainfall growing season could still occur, so lower 
yields are still possible. The probability distributions in Fig. 9 high¬ 
light the fact that under stochastic rainfall conditions, the question 
of optimization must be examined from a probabilistic point of 
view. Unlike in the previous section, optimal strategies for a sys¬ 
tem undergoing stochastic forcing must attempt to maximize profit 
while also being robust to adverse conditions, such as drought or 
flood years. This necessarily involves tradeoffs between maximiz¬ 
ing yield and profit on the one hand and mitigating risk on the other. 
For example, we point out that in the bottom of Fig. 9, the theoret¬ 
ical maximum (which refers to the value which would be obtained 



r[d] 


1.1 

1 

0.9 

0.8 

0.7 

0.6 

0.5 

0.4 


0 20 40 60 80 100 

r[d] 

Fig. 10. Crop yield as a function of £, the fraction of the total fertilization amount 
which is applied at the beginning of the growing season, and r, the time between 
the first and second fertilizer applications, for three soil depths: Z=33cm (top), 
Z= 67 cm (middle), and Z= 100 cm (bottom). In this figure, microirrigation was used 
to prevent the soil moisture from dropping below S*, and therefore the crop can 
experience only nitrogen stress, not water stress. 



if the crop experienced no water stress and took up as much N as 
possible) for profit (dashed line) occurs at a much lower fertiliza¬ 
tion rate than the maximum of the actual profit under stochastic 
rainfall conditions (black line), demonstrating in a simple way the 
necessity of accounting for the possibility of adverse conditions. 
A detailed analysis of such concerns is beyond the scope of this 
work, though we point out that many studies have examined the 
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related concept of resilience in ecological and social systems (see 
for example Walker et al. (2004)). 

In order to examine the impact of stochastic forcing on optimal 
fertilization, we first suppose that the total fertilization over the 
course of a growing season F t (the value of which can be found in 
Table 1 ) is to be divided into two treatments, which corresponds to 
a typical fertilization schedule for corn (e.g., Brady et al. (1996)). The 
placement of the fertilizer applications in time is varied by chang¬ 
ing the fraction of F t which is applied in each application and the 
amount of time between the two applications. In order to focus on 
optimal fertilization timing and amounts under stochastic condi¬ 
tions, we do not consider the other potential degrees of freedom 
in the fertilization scheduling, such as varying the total amount of 
fertilizer used or using more than two applications. The effect of 
varying £ and r (the fraction of F t in the first application and the 
time between the first and second applications, respectively) on 
crop yield can be seen in Fig. 10, which shows the yield response 
to varying £ and r for three different soil depths. Larger soil depths 
lead to less variation in 5 and thus less percolation and leaching, 
thereby increasing the fraction of £ - r space in which higher yields 
can occur. The exact location of the peak yield is a result of the bal¬ 
ance which maximizes the uptake of nitrogen and minimizes the 
loss due to leaching. 

6. Conclusion 

We have presented a dynamical system for crop evolution, based 
on the AquaCrop model (Steduto et al., 2009) and minimal mod¬ 
els for soil moisture and nitrogen cycling used in ecohydrology 
(Rodriguez-Iturbe and Porporato, 2004). It includes canopy cover, 
soil moisture, and soil nitrogen as its main state variables and tracks 
fluxes of water and nitrogen from evapotranspiration, nitrogen 
uptake, and leaching. This parsimonious model, with its reduced 
number of parameters, may be useful for evaluating the impact of 
different fertilization and irrigation strategies as well as different 
precipitation and climate regimes on crop yield, expected profit, 
and other outputs of interest. A simple objective function was used 
to compare optimal strategies of fertilization and irrigation. 

The results highlight the importance of considering, from a 
quantitative and theoretical point of view, the optimization of 
these agricultural inputs, and also provide a direct connection with 
climate parameters. Hydroclimatic forcing is a major driver of vari¬ 
ability in agricultural systems, which has implications not only for 
crop yield and profitability but also for environmental impact. The 
model developed here is capable of characterizing the variability 
in the model outputs and provides a link to the random processes 
which drive this variability. 

Agroecosystems cover a large portion of the Earth’s surface and 
provide essentially all of the global food supply. It is thus crucial to 
have a more complete understanding of the fluxes of water and 
nutrients in such systems, and their dependence on potentially 
changing hydroclimatic inputs and human activities. To this regard, 
the model presented here may be useful to explore scenarios and 
generate hypotheses. The framework can be extended in a number 
of directions. In order to emphasize the dynamical systems point of 
view, the model presented here necessarily included certain simpli¬ 
fications. However, including more detailed plant and soil models 
and performing a comparison with more complete crop models 
would provide firmer ground from which to make predictions. 
Moreover, the model could easily account for periodic seasonal 
variations in temperature, radiation, or rainfall, which alter the 
water and nutrient cycles and therefore the optimal fertilization 
and irrigation strategies. Finally, the nature of agroecosystems 
is that they are heavily intertwined with human activities (e.g., 
Sivapalan et al. (2012), Porporato et al. (2015), Assouline et al. 


(2015)), suggesting the need to couple models for ecological sys¬ 
tems and landscape evolution with social and behavioral models 
(e.g., harvesting in Parolari and Porporato (2016) and Pelak et al. 
(2016)). We hope that these considerations will be accounted for 
in future contributions, providing a quantitative framework for the 
sustainable use of soil and water resources while ensuring food 
security. 
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